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Abstract. We present the first numerical solutions of the 
coupled Einstein-Maxwell equations describing rapidly ro- 
tating neutron stars endowed with a magnetic field. These 
solutions are fully relativistic and self-consistent, all the ef- 
fects of the electromagnetic field on the star's equilibrium 
(Lorentz force, spacetime curvature generated by the elec- 
tromagnetic stress-energy) being taken into account. The 
magnetic field is axisymmetric and poloidal. Five dense 
matter equations of state are employed. The partial differ- 
ential equation system is integrated by means of a pseudo- 
spectral method. Various tests passed by the numerical 
code are presented. The effects of the magnetic field on 
neutron stars structure are then investigated, especially 
by comparing magnetized and non-magnetized configura- 
tions with the same baryon number. The deformation of 
the star induced by the magnetic field is important only 
for huge values of B (B > 10 10 T). The maximum mass 
as well as the maximum rotational velocity are found to 
increase with the magnetic field. The maximum allowable 
poloidal magnetic field is of the order of 10 14 T and is 
reached when the magnetic pressure is comparable to the 
fluid pressure at the centre of the star. For such values, 
the maximum mass of neutron stars is found to increase 
by 13 to 29% (depending upon the EOS) with respect to 
the maximum mass of non-magnetized stars. 
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1. Introduction 

Neutron stars are known to possess strong magnetic fields. 
Their polar field strength is deduced from the observed 



spin slowdown of pulsars via the magnetic dipole brak- 
ing model; for the 558 pulsars of the catalog by n Taylor 
et al. (1993), it ranges from B = 1.7 10~ 5 GT (PSR 
B1957+20) up to B = 2.1 GT (PSR B0154+61), with a 
median value B = 0.13 GT, most of young pulsars hav- 
ing a surface field in the range B ~ 0.1 — 2 GT. From the 
theoretical point of view, a considerable amount of studies 
have been devoted to the structure of the magnetic field 
outside the neutron star, in the so-called magneto sphere, 
in relation with the pulsar emission mechanism (for a re- 
view, see e.g. Michel 1991). The studies of the magnetic 
field inside neutron stars are far less abundant. Only re- 
cently, some works have been devoted to the origin and 
the evolution of the internal magnetic field, all of them in 
the non-relativistic approximation (Thompson & Duncan 
1993, Urpin & Ray 1994, Wiebicke & Geppert 1995, Urpin 
& Shalybkov 1995). 

Beside these studies of neutron star magnetic field, 
there exists a growing number of numerical computations 
of rapidly rotating neutron stars in the full framework of 
general relativity, taking into account the most sophisti- 
cated equations of state of dense matter to date (cf. Sal- 
gado et al. 1994a, b and reference therein, as well as Cook 
et al. 1994b, Eriguchi et al. 1994, Friedman & Ipser 1992). 
But in all these models the magnetic field is ignored. The 
present work is the first attempt to compute numerical 
models of rotating neutron stars with magnetic field in a 
self-consistent way, by solving the Einstein-Maxwell equa- 
tions describing stationary axisymmetric rotating objects 
with internal electric currents. In this way, the models 
presented below 

1. are fully relativistic, i.e. all the effects of general rela- 
tivity are taken into account, on the gravitational field 
as well as on the electromagnetic field. 
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1 In this article, we systematically use S.I. units, so that the 
magnetic field amplitude is measured in teslas (T) or, more 



conveniently in gigateslas (1 GT 
1 GT = 10 13 gauss. 



10 9 T). We recall that 
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2. are self-consistent, i.e. the electromagnetic field is gen- 
erated by some electric current distribution and the 
equilibrium of the matter is given by the balance be- 
tween the gravitational force, the pressure gradient and 
the Lorentz force corresponding to the electric current. 
Moreover, the electromagnetic energy density is taken 
into account in the source of the gravitational field. 

3. give the solution in all space, from the star's centre to 
infinity, without any approximation on the boundary 
conditions. 

4. use various equations of state proposed in the litera- 
ture for describing neutron star matter. 

The restrictions of our models are the following ones: 

1. We consider strictly stationary configurations. This ex- 
cludes magnetic dipole moment non-aligned with the 
rotation axis. Indeed, in the non-aligned case, the star 
radiates away electromagnetic waves as well as gravi- 
tational waves (due to the deviation from axisymmetry 
induced by a magnetic axis different from the rotation 
axis); hence it loses energy and angular momentum, so 
that this situation does not correspond to any station- 
ary solution of Einstein equation. Thus the stationary 
hypothesis implies that we restrict ourselves to axisym- 
metric configurations, with the magnetic axis aligned 
with the rotation axis. 

2. Moreover, we consider only poloidal magnetic fields 
(i.e. B lying in the meridional planes). Indeed, if the 
magnetic field had, in addition to the poloidal part, 
a toroidal component (i.e. a component perpendicular 
to the meridional planes), the circularity property of 
spacetime would be broken, which means that the two 
Killing vectors associated with the stationarity and the 
axisymmetry would no longer be orthogonal to a fam- 
ily of 2-surfaces (cf. Carter 1973, p. 159). In the circu- 
lar coordinate system (t, r, 9, 0) can be found 
such that the components of the metric tensor g are 
zero except for the diagonal terms and only one off- 
diagonal term (gt<f>)- In the non-circular case, only one 
component of g can be set to zero (g r d), resulting in 
much more complicated gravitational field equations 
(Gourgoulhon & Bonazzola 1993). On the contrary, 
perfect fluid stars with purely rotational motion (no 
convection) generate circular spacetimes and poloidal 
magnetic fields preserve this property (Carter 1973). 

3. Being not interested in modelling pulsar magneto- 
spheres, we suppose that the neutron star is sur- 
rounded by vacuumH. 

The numerical code we use is an electromagnetic exten- 
sion of the code presented in Bonazzola et al. 1993 (here- 
after BGSM), which was devoted to perfect fluid rotat- 
ing stars and used to compute neutron star models with 



2 by vacuum we mean that there is no matter outside the 
star; nevertheless there is some electromagnetic field, so that 
the total stress-energy tensor (right-hand side of the Einstein 
equation) is not zero outside the star (electrovac spacetime). 



various equations of state of dense matter (Salgado et 
al. 1994a, b). We will not give here the complete list of 
the equations to be solved but only the electromagnetic 
ones (Maxwell equations) in Sect. ||, referring to BGSM 
for the gravitational part. Likewise we will not present the 
numerical technique, based on a pseudo-spectral method, 
since it has been detailed in BGSM. We will discuss in 
Sect. H only the numerical procedure and tests of the elec- 
tromagnetic part of the code. We analyze the effects of the 
magnetic field on static configurations in Sect. |] and on 
rotating configurations in Sect. [| the Sect. ^ describing 
in detail the case of constant baryon number sequences. 
Finally Sect, fj] summarizes the main conclusions of this 
study. 

2. Electromagnetic equations 

2.1. Definitions and notations 

Following BGSM, we use MSQI {Maximal Slicing - Quasi- 
Isotropic) coordinates (t,r,9,(j>), in which the metric ten- 
sor g takes the form 

g a p dx a dx 13 = —N 2 dt 2 + A A \B 2 r 2 sin 2 9 (# - N^dtf 



■b,(* 



r 2 +r 2ja2\ 



2 d6 2 ) 



(1) 



where N, A and B are four functions of (r, 9), some- 
times represented by their logarithms: 



:=lniV 



:=h\A ; (3:=lnB. 



(2) 



Carter (1973, theorem 7) has shown that the most gen- 
eral form of the electric 4-current j compatible with the 
hypothesis of stationarity, axisymmetry and circularity 
has the following components with respect to the (t, r, 9, <fi) 
coordinates: j a — (j*, 0, 0, j^). He has also shown that the 
corresponding electromagnetic field tensor F must be de- 
rived from a potential 1-form A with the following com- 
ponents A a = (A t ,0, 0, A^): 



F a p — Ap :a — A a ^ 



(3) 



We will not use the "orthonormal" electromagnetic poten- 
tials $ and \I> introduced in BGSM but work directly with 
the components A t and A^. A nice feature of A^ is that 
the magnetic field lines lie on the surfaces A^ — const. 

For our problem, a priviledged observer is the observer 
Co whose 4-velocity n is the unit vector normal to the 
t = const hypersurfaces (hereafter referred to as the E t hy- 
persurfaces): it corresponds to the observer "at rest with 
respect to the star's centre" of the Newtonian physics. 
It is called the locally nonrotating observer by Bardeen 
(1970), the zero- angular-momentum observer (ZAMO) by 
Bardeen (1973), the Eulerian observer by Smarr & York 
(1978) and the FIDO by Thome et al. (1986). The elec- 
tric field E and the magnetic field B as measured by the 
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observer Oq are given b; 
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2.3. Lorentz force and condition for a stationary motion 
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(5) 



where t a f3iiv is the Levi-Civita tensor associated with 
the metric g. Note that a consequence of Eq. (5) is that 
the magnetic field lines lie on surfaces Ad, — const. In 
non-relativistic studies, Ad, is usually called the magnetic 
stream function or magnetic flux function and is denoted 
a or A (see e.g. Sakurai f985, Sauty & Tsinganos f994). 

2.2. Maxwell equations 

The source-free Maxwell equations -F] a /3 ;7 ] = arc au- 
tomatically satisfied by the form (|J) of F. The remain- 
ing Maxwell equations F a ^ = /io j a can be expressed in 
terms of A t and A$ as the Maxwell-Gauss equation 

A* , . t A 4 B 2 



A, A f 



-Mo-gj (9tt f + 9t4, j*) 



N 2 



-N't'r 2 sin 2 i 



xdAtdN 4, - 1 + 
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1 dA* 



dr r tan 6 d6 
and the Maxwell- Ampere equation 

A 3 I = -no A 8 (j 4 - ' 
A 4 B 2 



(6) 
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r sin i 



r sine dN^ (dA t 
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where 



A, 



(7) 



(8) 



da dp := 



r sin 9 

use has been made of the abridged notation 
dad£ \_dadj3_ 
dr dr r 2 dO d9 ' 

and A3 and A3 are, respectively, the scalar Laplacian and 
the 4> component of the vector Laplacian, in three dimen- 
sional flat space: 
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3 Unless otherwise specified, we will refer in the following to 
E and B as the electric field and the magnetic field, without 
mentioning "as measured by the observer Co" ■ 



A first integral of the equation of fluid stationary mo- 
tion (momentum conservation equation), taking into ac- 
count the Lorentz force exerted by the electromagnetic 
field on the conducting medium has been derived in BGSM 
(Eq. (5.30)); it reads 

H{r,0) +v(r, 6) - In T(r,6) + M(r, 6) = const. , (12) 
where 

— H is the fluid log-enthalpy defined as 
ff:=ln(^) , (13) 

e being the fluid proper energy density, p the pressure, 
n the proper baryon density, niB a mean baryon rest 
mass (mb = 1.66 10 -27 kg); 

— v is the gravitational potential defined by Eq. (||); 

— r is the Lorentz factor relating the observer Oo and 
the fluid comoving observer (cf. Eq. (5.6) of BGSM); 

— M is the electromagnetic term induced by the Lorentz 
force. 

M is actually a function of Ad, only, expressible as 



M(r,0)=M(Ad,(r,6)) = - 



f(x) dx 



(14) 







where / is an arbitrary regular function relating the com- 
ponents of the electric current to the electromagnetic po- 
tential Ad, as an integrability condition of the equation of 
fluid stationary motion (cf. Eq. (5.29) of BGSM): 



j*-ilj t =(e + p)f{A 4> ) , 



(15) 



where f2 is the fluid angular velocity defined as £1 :— 
vfi /vf. We will call / the current function. Different 
choices for / will lead to different magnetic field distri- 
butions. Note that in the limit of an incompressible New- 
tonian body, the integrability condition (|lj]) reduces to 
the relation derived by Ferraro (1954) [cf. his Eqs. (11) 
and (14), his function U being exactly our A^ and his 
function / being —p,Qpc 2 times our function /]. 



(9) 2.4- Perfect conductor relation 

According to Ohm's Law, and assuming that the matter 
has an infinite conductivity, the electric field as measured 
by the fluid comoving observer, E' a = F a pu^, must be 
zero. This condition leads to the following relation be- 
tween the two components of the potential 4- vector inside 
the star: 



dA t 

dx 1 



-n 



dA d 
dx 1 



(16) 



As shown in BGSM, a stationary configuration with a 
magnetic field is necessarily rigidly rotating (i.e. has £1 = 



3 
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const.). Consequently, Eq. jl^ ) is integrated immediately 
to 



A t = -£IA 6 + C 



(17) 



where C is a constant. The choice of C will fix the total 
electric charge of the star. Note that Eq. ( |l~7j ) holds inside 
the star only. 

2.5. Global quantities 

The magnetic dipole moment M. is given by the leading 
term of the asymptotic behaviour of the magnetic field as 
measured by 0$: 



(r) — >oc 



fio M cos I 
2n r 3 



B 



(0) ^r- *oo 



fio A4 sin 9 
47T r 3 



. (18) 



Bi r ) and Brg\ being the component of B in the orthonor- 
mal basis associated to (r,6,(f>). They are related to the 
components B r and Bg given in Eq. (5) by 



B 



B 



(19) 



The total electric charge of the star, Q, is given by the 
leading term of the asymptotic behaviour of the electric 
field as measured by Oq: 



E, 



(r) 



fiO_Q_ 

47T r 2 



E, 



(S) 



Equivalcntly, Q is expressible as 



A*o Q_ 

47r r 2 



NA e 

— — j* r 2 sin 9 dr alO d<j> + Q S urf , 
B 

where Q SU ri is the surface charge. 



(20) 

(21) 
(22) 



3. The numerical procedure and its tests 

We describe here only the electromagnetic part of the 
code. For the fluid and gravitational part, we refer to 
BGSM. Let us simply recall that our code is based on an 
iterative procedure, each step consisting in solving Pois- 
son equations by means of a Chebyshev-Legendre spec- 
tral method developed by Bonazzola & Marck (1990). 
Typically the computations have been performed with 41 
Chebyshev coefficients in r and 21 coefficients in 9. 

3.1. Numerical resolution of the electromagnetic equations 

In order to obtain a solution of the equations presented 
in Sect. @, one needs to specify a value for the central 
log-enthalpy H, a value for the angular velocity Q, and a 
value for the total electric charge Q (possibly Q = 0). (An 
alternative scheme could be, instead of fixing Q, to impose 
EB = as a boundary condition at the surface of the star 



in order to mimic a magnetosphere; the solution would be 
then an approximate one outside the star, which is not the 
case in the procedure described below) . One needs also to 
pick a choice for the current function /. 

The iterative procedure is then as follows. At the first 
step, A^ is set identically to zero. At a given step, a value 
°A t of A t inside the star is deduced from the value of A^ 
at the previous step by the perfect conductivity equation 
(p7|), setting the constant C to zero. Then one computes 
the Laplacian of °A t and use the Maxwell-Gauss equation 
(6) to obtain the value of the charge density j l inside the 
star (using for the which appears in Eq. (6) its value 
at the previous step). The current density component 
is then deduced from j* and the previous value of A^ via 
the function / and Eq. ( |l5| ) . Then one solves the Maxwell 
equations (6) and (7) by considering them as two linear 
Poisson equations for A t and A^, with a fixed right-hand 
side: this latter is calculated from the newly determined 
f and and the previous step values of A t and A^ (in 
the non- linear terms containing dA t and dA^). For A^, 
the boundary conditionQ A^ — at r = +oo is used. We 
then obtain a unique smooth solution for A$. For A t the 
procedure is a little more complicated since in general a ro- 
tating perfect conductor is endowed with a surface charge 
density, so that the component of the electric field normal 
to the surface is discontinuous, meaning that the electric 
potential A t is not differentiable accross the surface. We 
then solve the equation for A t in two steps: (i) a solution 
1 A t to the Poisson equation (6) is obtained outside the 
star with the boundary condition At — at r — +oo. 
This solution, is by no means unique, since any harmonic 
function ip (i.e. satisfying A^ip = 0) vanishing at infinity 
can be added to At- (ii) an harmonic function of the form 



ii a\ ST P i( cost 



(23) 



1=0 



where Pi is the Legendre polynomial of degree I, is added 
to x At ouside the star so that °A t := x A t + ip matches, on 
the star's surface, the value of °A t inside the star which has 
been determined previously by the perfect conductor re- 
lation °At = — flA<p. Technically the matching is achieved 
by solving a linear system for the L + 1 coefficients ai of 
the expansion (^3|) , L + 1 being the number of grid points 
on the star's surface (number of collocation points in 9). 
The solution °A t obtained in this way is continuous in all 
space but not necessarily differentiable accross the star's 
surface. Besides, this solution has a certain electric charge, 
°Q say, given by Eq. ( pfj| ) and which does not a priori co- 
incide with the desired electric charge Q. The arbitrary 
nature of the constant C of Eq. (|i~7| ) can be used to adjust 
the electric charge, °At having been obtained by setting 



4 as explained in BGSM, our numerical grid extends to r = 
+oo thank to a change of coordinate u = 1/r in the space 
outside the star. 
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C = 0. This is achieved by considering the following so- 
lution of the equation A3 2 At — 0: 2 At — 1 inside the star 
and 2 At has the same form as the function ifi of Eq. ( p3| ) 
outside the star, the coefficients ai being determined in 
order to insure the continuity of 2 A t accross the star's sur- 
face. Let be the electric charge associated with 2 At. 
The constant C is then determined by 

and the final solution for A t is 

A t = °A t + C = -SIA^ + C inside the star (25) 
A t = °A t + C 2 A t outside the star . (26) 

By construction A t satisfies the Poisson equation (6) (with 
the right-hand computed from values at the previous step) 
and corresponds to the required electric charge Q. 

Once At and A^ are obtained, the electric and mag- 
netic fields are computed via Eqs. (4) and (5). The elec- 
tromagnetic stress-energy tensor is then deduced from 
Eq. (5.32) of BGSM and put as a source term, beside the 
fluid stress-energy tensor, in the Einstein equations giving 
the gravitational field (Eqs. (3.19)-(3.22) of BGSM). 

3. 2. Tests of the numerical code 

Various tests have been performed on the electromagnetic 
part of the code (for the gravitational part we refer to 
BGSM), mainly by comparing numerical solutions to ex- 
act ones for simple charge and current distributions. 

3.2.1. Rotating magnetic dipole 

We have considered an electrostatic dipole, a magneto- 
static dipole and a rotating magnetic dipole, all of them 
in flat space (no gravitation). Let us detail the results for 
the rotating magnetic dipole, since they are representative 
of the precision achieved in all other cases. 

By "rotating magnetic dipole" , we mean two differ- 
ent physical situations: (i) a conducting sphere in rotation 
with surface electric currents and (ii) a conducting sphere 
in rotation with an infinitely small electric current loop at 
the origin. In case (i) , the current distribution is modelled 

by 

jffl (r, 6) = jb v sin 9 exp (- (? ~ ) , (27) 

where jo is a constant, R is the sphere's radius and 6 is 
the width of the distribution, which we do not put to zero 
in order to have only smooth functions. Typically, we use 
5 ~ .01 R. The exact solution for (E, B) (corresponding 
to 5 = 0) may be found in Ruffmi & Treves (1973). The 
magnetic field is uniform and parallel to the rotation axis 
inside the sphere and dipolar outside. The electric field 



5 

outside the sphere is the sum of a monopolar term and a 
quadrupolar one. 

In case (ii), the current distribution is modelled by 

j W (r,9)=j rsm6eiq ? (- 1 ^J . (28) 

The exact solution (corresponding to S = 0, i.e. to a point 
magnetic dipole moment) may be found in the Table 4.1 
of Michel (1991)El. The magnetic field is dipolar in all 
space. The electric field outside the sphere is the sum of 
a monopolar term and a quadrupolar one. The difference 
between the exact solution and the solution given by the 
numerical code is depicted in Fig. ^ in the case of a con- 
ducting sphere of diameter 1 m, rotating at the angular 
velocity f2 = 3000 rads -1 , with a central loop current 
jo = 10 11 Am" 2 . Near the star's centre, the discrepancy 
is quite important because the current distribution used 
in the code is not exactly a point dipole (S ^ 0). But at 
half the stellar radius, the relative error becomes less than 
10 -5 and is of the order ~ 10 -9 outside the star. 



Relative error 01 


electric field 




Relative error on magnetic field 

















CD 













5x10 4 10 3 1.5x10 3 5x10 4 10 3 1.5x10 3 



Coordinate radius [km] Coordinate radius [km] 

Fig. 1. Comparison between the exact solution and the numer- 
ical one in the case of a rotating conducting sphere with a point 
magnetic dipole moment at its centre. The left figure corre- 
sponds to the electric field, the error being defined as the square 
root of [(E (r) - Effi*)* + (E (e) - E^f]/[{E {r) f + (E m ) 2 ]. 
The right figure depicts the same thing for the magnetic field. 
Asterisks denote the sphere's surface. 



3.2.2. Comparison with Ferraro's solution 

Ferraro (1954) has obtained an analytical solution for a 
Newtonian incompressible fluid body with an electric cur- 
rent distribution corresponding to the choice f(x) — const 

5 in this reference, the term 2/3r in the seventh line should 
read 2/3r 2 
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for the current function of Eq. (|l5|). Ferraro's solution is 
valid for a body whose shape remains close to a sphere 
and reads 

A$ = -fi pc 2 fo(r 2 /l0 - R 2 /6)r 2 sin 2 9 for r<R (29) 



=tt/2 — , 0=V4 -~ 0=0 



R 5 

= ^pc 2 f Q —sm 2 ( 
15r 



for r > R 



(30) 



where /o is the constant value taken by the current func- 
tion / and R is the stellar radius. The discrepancy between 
Ferraro's solution and the numerical one resulting from 
the code is depicted in Fig. |[ It is everywhere less than 
10 -3 , being less than 5 x inside the star. These num- 



bers are higher than those of the test of Sect. |3.2.1| but one 
should keep in mind that Ferraro's solution is only an ap- 
proximate one, assuming no deviation from spherical sym- 
metry, whereas the numerical solution takes into account 
the deformation of the star by Lorentz forces. Of course, in 
order to perform the test, we have considered a magnetic 
field weak enough for Ferraro's approximation to be valid, 
as well as a very weak gravitational field for the Newtonian 
approximation to be justified: the tested configuration is a 
2.67xlCT 2 M star of (constant) density 1.66 xlO 9 kgm~ 3 
with a polar magnetic field i? po i c = 6.6 x 10 5 T. 

Ferraro's solution is interesting since it allows to test 
the response of the star to the magnetic field, contrary 
to the tests of Sect. 3.2.1. The deformation of the star as 
computed by Ferraro at first order around the spherical 
symmetry, is given in terms of the polar magnetic field by 
the eccentricity 



15 B 



pole 



4 y/irGfio pR 



(31) 



For the configuration considered above, e = 0.04670. The 
eccentricity resulting from the code is 0.04683, so that 
the relative discrepancy is 3 x 10 -3 . Here again let us re- 
call that Ferraro's solution is only an approximate one, so 
that this value reflects not only the true numerical error 
but also Ferraro's error with respect to the exact solution. 
Nevertheless the very good agreement with Ferraro's so- 
lution means that the action of Lorentz forces upon the 
fluid is correctly treated by the code and in particular that 
there is no error in the conversion from magnetic units to 
fluid units. 

3.2.3. Virial identities 

A different kind of test is provided by the virial identi- 
ties GRV2 (Bonazzola & Gourgoulhon 1994) and GRV3 
(Gourgoulhon & Bonazzola 1994), the latter being a rel- 
ativistic generalization of the Newtonian virial theorem. 
GRV2 leads to the |1 — A| error indicator introduced in 
BGSM. These virial identities allow to test each compu- 
tation and not only the simplified ones presented above. 
It notably controls the convergence of the iterative proce- 
dure used in the code. In most of our calculations, the final 




Fig. 2. Comparison between Ferraro's analytical solution and 
the numerical one in the case of a Newtonian incompressible 
fluid endowed with a magnetic field corresponding to the cur- 
rent function f(x) — const. The plotted quantity is the relative 
difference between the two values of A$ as a function of the 
radial coordinate r for three values of 8. Asterisks denote the 
star's surface. 



value of the GRV2 or GRV3 error indicator is of the order 
10~ 5 (polytropic EOS) or a few 10" 4 (tabulated EOS). 
Let us recall that tabulated EOS introduce an important 
numerical error because they do not stricly obey to the 
thermodynamical relations (cf. Sect. 4.2 of Salgado et al. 
1994). We systematically rejected any solution for which 
the GRV2 or GRV3 error is greater than 10 -2 , such a 
relative error being too large for considering that the nu- 
merical procedure has converged. This may happen if the 
magnetic field is too large for any stationary configuration 
to exist. 



4. Static magnetized configurations 

In order to investigate some purely magnetic effects on 
stellar configurations, let us first consider the case of static 
(i.e. non rotating) neutron stars. More precisely, a static 
spacetime (in the sense of the time Killing vector being 
hypersurface-orthogonal) implies not only a non rotating 
fluid (Q = 0) but also a vanishing electric charge. Only in 
this case the electric field as measured by Oo is zero. On 
the contrary, a non-vanishing electric charge creates an 
electrostatic field outside the star, so that there exists a 
net electromagnetic momentum density vector E x B/2/irj 
which generates a non-zero angular momentum, though 
the fluid does not rotate [Feynman's disk paradox (Feyn- 
man et al. 1964, Ma 1986, de Castro 1991)]. Consequently 
the shift vector component ^ and the spacetime is 
not static. We observed this effect numerically but only 
for huge values of the electric charge: Q ~ 10 17 C. 
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-20 20 

x [ km ] 

Fig. 3. Fluid proper density isocontours in the (r, 9) plane of 
a static magnetized star of 2.98 Mq built on the Pol2 EOS of 
Salgado et al. (1994a) (polytropic EOS with 7 = 2) with a 
current function f(x) = const. = 4 x 10 15 /i? Am _2 p~ u 1 c c~ 2 , 
where R is the r coordinate of the star's equator. The central 
energy density is e = 1.42 p nuc c 2 (1 p nuc := 1.66 x 10 17 kgm -3 ). 

4-1- Deformation of the star 

Even in the static case, the stress-energy tensor of the elec- 
tromagnetic field is not isotropic (cf. Eq. (5.33) of BGSM 
with U = 0) so that the star deviates from spherical sym- 
metry. For large magnetic field this deviation becomes im- 
portant as shown in Fig. [| The corresponding magnetic 
field lines are plotted in Fig. |]. Note that the situation de- 
picted in Figs. |^ and ^| corresponds to an important mag- 
netic field: .Bpoic = 9.1 x 10 3 GT. The flattened shape is 
a general result for a magnetized configuration. We thus 
may conclude that the Lorentz forces exerted by the elec- 
tromagnetic field on the conducting fluid behave as cen- 
trifugal forces. Note that the flattening property of the 
magnetic field has been recognized long ago (in the New- 
tonian case) by Chandrasekhar & Fermi (1953). 

4-2. Maximum values of the magnetic field 

If the value of the central magnetic field is high enough for 
the magnetic pressure equals the fluid pressure p, then the 
total stress tensor Sij has a vanishing component along the 
symmetry axis, as it can be seen by combining Eq. (5.7c) 
and (5.32c) of BGSM: 

S r r =p+^-{B e B e -B r B r ) , (32) 

with, on the symmetry axis near the centre, BgB e = and 
B r B r /2fi ~ p. Hence S r r = 0. Away from the centre, 
the fluid pressure decreases more rapidly than the mag- 
netic pressure along the symmetry axis, so that S r r < 0, 




-50 

x [ km ] 



Fig. 4. Magnetic field lines in the (r, 8) plane for the 
static star configuration corresponding to Fig. ^. The thick 
line denotes the star's surface. The magnetic field ampli- 
tude is B c = 3.57 x 10 4 GT at the star's centre and 
-Bp i e = 9.1 x 10 3 GT at the north pole; the magnetic dipole 
moment is M = 5.42 x 10 32 Am 2 . 



o D" 
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-20 -10 10 20 

x [ km ] 

Fig. 5. Fluid proper density isocontours in the (r, 9) 
plane for the maximum mass static magnetized star 
built on the Pol2 EOS of Salgado et al. (1994a) 
(polytropic EOS with 7 = 2). The current function 
/(as) = const. = 5.63 x 10 15 /ii Am _2 p~ u 1 c c~ 2 , where R is the 
r coordinate of the star's equator. The central energy density 
is e — 1.42 pnucC 2 and the mass M = 4.06 Mq; the other char- 
acteristics are listed on the third line of Table pi. 
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Fig. 6. Magnetic field lines in the (r, 6) plane for the 
maximum mass configuration corresponding to Fig. |j| The 
thick line denotes the star's surface. The magnetic field am- 
plitude is B c — 9.00 x 10 4 GT at the star's centre and 
Bpoie = 4.57 x 10 4 GT at the north pole; the magnetic dipole 
moment is M = 1.12 X 10 33 Am 2 . 



Fig. 7. Ratio of the magnetic pressure B 2 /2/io to the fluid 
pressure p as a function of the radial coordinate r and along 
three angular directions for the maximum mass configuration 
displayed in Figs. ^ and ^. Outside the star, where p — 0, the 
logarithm of the ratio has been set to 5. At the star's centre, 
(B 2 /2fi )/p = 1.36. 



which means that the fluid + magnetic field medium devel- 
ops some tension instead of pressure along the symmetry 
axis. The star displays then a pinch along the symmetry 
axis, as it can be seen on Figs. || and ^. The fact that in 
this case the magnetic pressure exceeds the fluid pressure 
everywhere on the symmetry axis can be clearly seen in 
Fig. ^ (dotted line), which represents the ratio of the two 
pressures throughout the star. For magnetic fields higher 
than that presented in Figs. Jjj-jj], no stationary configura- 
tion can exist and the numerical procedure described in 



Sect. 3.1 fails to converge. 



The stellar shape displayed in Fig. || is a typical fea- 
ture of the maximum allowable magnetic field and is not 
due to the fact that this particular configuration is the 
maximum mass one or that the EOS is Pol2. Indeed, we 
found very similar shapes for all central densities in the 
explored range (0.15/? nuc up to 5p nuc ) and all EOS when 
the magnetic field has reached the maximum value for 
which a solution could be obtained. In all cases, the ratio 
of B 2 /2fio to p was close to 1 at the star's centre (cf. the 
sixth column of Table |l|) . 

Note that the maximum values of i? po ic are four orders 
of magnitude higher than the maximum observed value: 
B polc = 2.1 GT in the pulsar PSR B0154+61 (Taylor et al. 
1993), which means that the magnetic field in the observed 
pulsar is not limited by the condition of stellar equilibrium 
but by the physical processes it originates from. 



4-3. Effect of the magnetic field on M max 

In order to investigate the effect of the magnetic field on 
the maximum mass of static configurations, we choose the 
simplest form of the current function: 



f(x) = const. = /o , 



(33) 



and vary fo from (no magnetic field) up to the max- 
imal value for which the code converges and for which 
the ratio of the magnetic pressure to the fluid pressure at 
the star's centre come close to one (cf. sixth column of 
Table |l|) . The choice (^3|) results in electric currents con- 
centrated deep inside the star, with the same behaviour 
as that represented in Fig. [l(] for the rotating case. Each 
static configuration is determined by the value of the cen- 
tral log-enthalpy H c and the value of fo. In the plane 
(iJ c ,/o), we followed curves of constant magnetic dipole 
moment M and determined for each of them the maximal 
gravitational mass. 

This study has been conducted for five of the EOS used 
in Salgado et al. (1994a,b): Pol2 (7 = 2 polytrope), Di- 
azll (model II of Diaz Alonso 1985), PandN (pure neutron 
model of Pandharipande 1971), BJI (model IH of Bethe 
& Johnson 1974) and HKP (Haensel et al. 1981). For de- 
tails about these EOS, the reader may consult Sect. 4.1 of 
Salgado et al. (1994a). 

Some results are presented in Table [j] which gives max- 
imum mass configurations (i) with no magnetic field (first 
line), (ii) with some fixed magnetic dipole moment (sec- 
ond line), (iii) among all magnetized configurations (third 
line). 
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Different behaviours appear for different EOS: the 
maximum gravitational mass at fixed magnetic dipole mo- 
ment M as well as the maximum baryon number at fixed 
M. are both increasing functions of M. for the EOS Di- 
azll (cf. Fig. ^]), whereas the maximum baryon number is 
a decreasing function of M for the EOS BJI (cf. Fig. |9|), 
the maximum mass remaining a increasing function. Pol2 
models behave as DiazII ones, whereas PandN and HKP 
models behave as BJI ones. In all cases, the maximum 
gravitational mass increases with the magnetic field, by 
28.6% for Pol2, 20.2% for DiazII, 14.9% for PandN, 17.3% 
for BJI, and 13.3% for HKP (compare first and third lines 
in Table 0). 

5. Rotating magnetized configurations 

5.1. Electric charge 



fl=rr/2 



9=V 4 



= 



As we have seen in Sect. || and 3d, the total electric charge 
Q is a freely specifiable parameter of our models. Realistic 
neutron stars certainly possess a net electric charge, which 
may be of the order of 10 12 C (Michel 1991, Chap. 4 and 
Cohen et al. 1975). The key point is that this (huge) charge 
is compensated by the charge provided by the magneto- 
spheric particles to lead to globally neutral star + mag- 
netosphere system. Now, as stated in Sect, [l], we do not 
attempt to model any filled magnetosphere and consider 
neutron stars surrounded by matter vacuum. In conse- 
quence, we do not try to give a specific charge to our mod- 
els (this would require a detailed magnetosphere model) 
and consider globally neutral stellar configurations: Q = 0. 
Anyway, the effect of a non-zero value of Q on the non- 
electromagnetic characteristics of the star reveals to be 
extremely small: we verified that 



V H c ,n, \M(H C ,Q,Q = 10 18 C) - 
M(H c ,Q,Q = 0)\ < 



1(T 5 M m 



(34) 



which is of the order of the numerical error committed by 
the code. The global quantity which is the most sensitive 
to a non-zero value of Q is the angular momentum as 
mentioned in Sect. [|. In the rotating case however the fluid 
angular momentum dominates by at least one order of 
magnitude the electromagnetic angular momentum, even 
for huge values of Q (Q ~ 10 18 C). 

5.2. Structure of the electromagnetic field 

In order to consider different electric current distributions, 
concentrated around the star's centre or not, we use dif- 
ferent choices for the current function /. The most imme- 
diate choice is ( |33|) (f{x) = fo = const). The resulting 
electric current distribution is shown in Fig. [l^ for a Pol2 



M = 3.37 M Q model rotating at ft = 3 x KT rads" 



has a maximum at one third of the stellar radius. The 
resulting magnetic field is shown in Fig. 11 and O. 



The 

magnetic field amplitude has a maximum near the star's 




10 20 30 

Coordinate radius [ km ] 



Fig. 10. Electric current distribution induced by the current 
function f(x) = const. = 10 15 /i? A m _2 p~ u 1 c c~ 2 for aPol2 EOS 
model, of central energy density e c = 3.06 p mc c 2 , rotating at 
Q = 3 x 10 3 rads" 1 . The plotted quantity is the value j (<w of 
the azimuthal component of j in an orthonormal basis in the 
equatorial plane 8 = n/2 (solide line), along = 7r/4 (dashed 
line) and along the rotation axis 9 — (dotted line). Asterisks 
denote the star's equatorial radius. 



centre, the surface value being approximatively four times 
lower. The dominant dipolar structure of B clearly ap- 
pears in Fig. [l3|, as if B was generated by a single current 
loop located inside the star. The induced electric field is 
represented in Figs. [l2| and [b|. Note on Fig. |l2| the discon- 
tinuity in the component of E normal to the star's surface, 
due to surface charges. The cusp across the surface in the 
component Eiq\ along the radius 9 = 7r/4 is due to the dis- 
crepancy between Erg\ and the component of E tangential 
to the surface (resulting from the non-spherical shape of 
the star), this latter being perfectly smooth, as expected. 
The fact that A t is continuous but not differentiable across 
the star's surface, because of the surface charges, is clear 
on Fig. |l3|. Note also on Fig. |l3| that, inside the star, the 
lines A t = const coincide with the magnetic field lines, as 
expected from the perfect conductor relation ( |l7| ) and the 
fact that the magnetic field lines are the lines = const, 
(cf. Sect.|2~l|). 

In order to get current distributions different from the 
one described above, we have considered the following 
forms of the current function, beside the constant form 



f{x) 



l + x 
a(l- 



A° 



[fix) 



A^_ 

Al ' 



(35) 
(36) 
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Table 1. Maximum mass (at fixed magnetic dipole moment) neutron stars in the non-rotating case. For each EOS, the last line 
is the maximum mass configuration among all static magnetized models. M is the magnetic dipole moment, B c and £> po i c are 
respectively the values of the magnetic field at the centre and the north pole of the star, H c is the central value of the log-enthalpy 
as defined by Eq. (|l|), e c the central energy density (1 p nuc := 1.66 10 17 kgm 3 ), p m g,c/pfl,c the ratio of the magnetic pressure 
to the fluid pressure at the star centre, M the gravitational mass, B the baryon mass, J5bind the binding energy per baryon 
(with the convention -Ebmd < for a bound configuration), -Rcirc the circumferential radius (length of the star's equator divided 
by 2n) and GRV2 and GRV3 the estimates of the global numerical relative error provided by the virial identities GRV2 and 
GRV3 (cf. text). We use G = 6.6726 10~ n m 3 kg _1 s- 2 and 1 M© = 1.989 10 30 kg. 



EOS 
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-Be 


^polc 


H c 


e c 
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M 


B 
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-Rcirc 
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Am 2 ] 


[10 3 GT] 


[10 3 GT] 
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Pol2 











0.491 


4.17 





3.158 


3.470 


-0.0897 


21.77 


1E-14 


4E-12 




2.00 


37.8 


6.2 


0.483 


4.07 


0.04 


3.182 


3.486 


-0.0872 


21.93 


1E-06 


2E-06 




11.22 


90.0 


45.7 


0.225 


1.42 


1.36 


4.062 


4.279 


-0.0507 


26.45 


1E-03 


7E-04 


DiazII 











0.610 


15.21 





1.933 


2.210 


-0.1253 


10.92 


1E-04 


1E-04 




0.50 


73.8 


14.1 


0.600 


14.85 


0.03 


1.940 


2.212 


-0.1229 


11.00 


2E-04 


2E-04 




3.08 


207.4 


108.3 


0.285 


6.22 


1.10 


2.324 


2.508 


-0.0734 


13.12 


1E-04 


9E-05 


PandN 











0.733 


24.39 





1.662 


1.932 


-0.1340 


8.55 


1E-04 


2E-04 




0.20 


64.7 


11.9 


0.727 


24.13 


0.01 


1.663 


1.931 


-0.1384 


8.57 


1E-04 


2E-04 




1.86 


302.8 


153.8 


0.350 


11.21 


1.01 


1.910 


2.064 


-0.0746 


10.00 


3E-04 


2E-04 


BJI 











0.699 


18.64 





1.856 


2.134 


-0.1303 


9.91 


2E-06 


3E-06 




0.30 


61.5 


11.0 


0.692 


18.40 


0.01 


1.858 


2.132 


-0.1288 


9.94 


2E-06 


2E-06 




2.63 


233.2 


121.0 


0.300 


7.47 


1.10 


2.176 


2.344 


-0.0717 


12.05 


1E-04 


5E-05 


HKP 











0.725 


8.75 





2.836 


3.422 


-0.1712 


13.67 


1E-04 


8E-05 




0.80 


56.0 


12.4 


0.714 


8.57 


0.02 


2.840 


3.417 


-0.1689 


13.75 


1E-04 


9E-05 




4.90 


196.6 


104.5 


0.370 


4.36 


0.90 


3.212 


3.594 


-0.1063 


15.85 


6E-04 


3E-04 
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Fig. 8. Maximum gravitational mass (left) and maximum baryon number (right) along static sequences at fixed magnetic dipole 
moment as a function of the magnetic dipole moment for the DiazII EOS . 
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n = o, f =c 



[10 30 A.m Z ] 




nt [10 30 A.m 2 ] 



Fig. 9. Same as Fig. g but for the BJI EOS 





Fig. 13. Magnetic field lines (left) and electric isopotential lines At = const (right) in the (r, 9) plane for the configuration 



Oef, 

considered in Figs, hd, [ill and Il2[ The thick line denotes the star's surface. 



where A , = Air x 10 4 i? 3 Tm , R being the value of a cyclic and/or chaotic behaviour when / is non- linear and 



the coordinate r at the equator. One might think about 
other choices, such as f(x) — ax + (3, f(x) — ax 2 , or 
f(x) = a/(l + x 2 ), but with such functions, the iterative 



procedure described in Sect. 3.1 revealed not to converge, 
resulting in a solution either cyclic or chaotic with re- 
spect to the number of steps. This behaviour is not due to 
some numerical instability but results from the non-linear 
coupling between and A^ via the Eqs. (7) and (|l~5|), 
which may be written schematically as a unique equation 
AAfj, — F(A ( j > ). Now it is well known that some sequences 
defined by an iterative relation u n+ i = F(u n ) may exhibit 



not one-to-one [take F(x) — Aa{x — x 2 ), which generates 
cycles of increasing order for 0.75 < a < 0.892 and chaotic 
sequences for a > 0.892]. For this reason, we did not con- 
sider the current functions mentionned above in our study. 



The choice (35) leads to a current distribution slightly 
different from that corresponding of f(x) — const. 
(Fig. 0), being simply more concentrated towards the 
star's centre, especially in the direction of the equatorial 
plane. The resulting electromagnetic field is very similar 
to that of Figs. Oy. 
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Coordinate radius [ km ] Coordinate radius [ km ] Coordinate radius [ km ] 



Fig. 11. Components _B( r j (left) and Bjg) (right) of the mag- 
netic field generated by the electric current distribution of 
Fig. [j^. The notations are the same as in Fig. [l^. 



e-n/2 — , e=n/4 - 




10 ao 30 

Coordinate radius [km] 



8=tt/2 — , 0=n/i - 



10 20 30 

Coordinate radius [km] 



Fig. 12. Components E^ (left) and (right) of the electric 
field corresponding to the configuration of Figs. [To] and [Tl| The 
notations are the same as in Fig. On the right figure, the 
dotted line merges with the solid one. 



Fig. 14. Electric current distribution induced by the choice 
f(x) = 10 4 [1 - l/((3.3 x 1(T 3 xf + 1)]/R Am^VnucC -2 for a 
Pol2 EOS model, of central energy density e c = 3.06 p nU cC 2 , 
rotating at SI = 3 x 10 3 rads -1 . The notations are the same as 
in Fig. nfj. 



8=7r/2 — , S=tt/4 - 
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Coordinate radius [ km 



9=tt/2 — , 6=7r/4 - 
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Fig. 15. Components Br r ) (left) and Big) (right) of the mag- 
netic field generated by the electric current distribution of 
Fig. [H| The notations are the same as in Fig. [h| 



The choice (36) leads to a current distribution which 
is represented in Fig. [l4|. There is no current around the 
centre and the maximum of the distribution is reached at 
almost half of the stellar radius. The resulting magnetic 
field is shown in Fig. [l^. The difference with that gener- 
ated by f(x) = const (cf. Fig. |ll|) is not important; sim- 
ply the maximum of B is achieved apart from the star's 
centre (except near the rotation axis). The corresponding 
electric field is shown in Fig. Ha It has more or less the 
same structure as that of Fig.ll2. 



5.3. Differences between magnetized and non-magnetized 
configurations 



For a given EOS and a given choice of the current function, 
we have explored the plane (H c , Q) by computing around 
14 x 14 = 196 models, with and without electromagnetic 
field. We present in this section the results concerning the 
effect of the electromagnetic field in the (H c , fl) space for 
the particular case of the Pol2 EOS and the form (33) 
of the current function, the constant /o being adjusted in 
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s=tt/2 — , 8=jt/4 - e=o ■ 



e=7r/2 — , 8=tt/4 — , 6=0 ■ 
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Coordinate radius [ km ] 
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Fig. 16. Components Er r ) (left) and Erg) (right) of the electric 
field corresponding to the configuration of Figs. ^ and [L5| The 
notations are the same as in Fig. [H]. On the right-hand figure, 
the dotted line is merged with the solid one. 




0.008 
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0.34 
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0.38 
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0. 



0.5 



Fig. 17. Difference between the baryonic mass [in M@] of con- 
figurations with a magnetic dipole moment M = 1.5 10 32 Am 2 
and the baryonic mass of configurations without magnetic field, 
as a function of the central log-enthalpy H c and the angular 
velocity Q [in 10 3 rads -1 ], for the Pol2 EOS. 



each model in order to lead to the same magnetic dipole 
moment M. — 1.5 10 32 Am 2 , so that the magnetized con- 
figurations may be considered as equivalent from a mag- 
netic point of view for an observer sufficiently far from the 
star. 

A first effect of the magnetic field is shown on Fig. [l?]: 
it increases the baryonic number at fixed H c and Q. This 
means that Lorentz forces, thanks to their centrifugal ef- 
fect (cf. Sect. 4.1), help the star to support more baryons 




Fig. 18. Same as Fig. [It] but for the total angular momentum 
J [in Mq c km] . 




Fig. 19. Same as Fig. [l?] but for the equatorial circumferential 
radius _R c irc fin kml. 



at the same central density and same rotation rate. It ap- 
pears on Fig. [I?] that this effect is rather insensitive to il, 
reflecting that this is the magnetic field which contributes 
the most to the Lorentz force, not the electric field. It can 
be verified that, as in the static case, the centrifugal effect 
of Lorentz forces enhances the oblateness of the star. 

A second effect of the magnetic field is to increase the 
total angular momentum J, as shown in Fig. [l^, more 
or less linearly with respect to f2. This augmentation is 
rather significative and can reach more than 1.5% of the 
zero magnetic field value. The redistribution of baryons 
inside the star induced by Lorentz forces is responsible for 



13 



14 



M. Bocquet et al.: Rotating neutron star models with magnetic field 




Fig. 20. Same as Fig. [l?] but for the equatorial fluid velocity 
as measured by the locally non-rotating observer, U aq [in units 
ofc]. 



most of the increase, well before the contribution to J of 
the electromagnetic field intrinsic angular momentum. 

The variation of the star's equatorial radius i? c i rc (de- 
fined as the star's equatorial circumference as measured 
by the observer Oq divided by 2ir) is depicted in Fig. |l9|. 
For small angular velocities, the radius increases, in the 
continuation of the static case where it has been seen 
that Lorentz forces stretch the star out (Sect. 4.1). But 



for larger angular velocities the radius decreases instead. 
This may be explained in the following way. At large rota- 
tion rates the stellar equilibrium is governed by the inertial 
centrifugal forces. Now these forces increase with the ra- 
dius and are rather important at the periphery of the star. 
On the contrary, the Lorentz forces are important near the 
centre and support there the massive nucleus of the star, 
so that for the same central density and rotation rate, the 
star can have a smaller radius to maintain the equilibrium, 
which is energetically favorable in so far as the peripheric 
inertial forces are near the shedding limit, contrary to the 
low rotation velocity case. Due to this reduction of the 
star's radius, the equatorial fluid velocity as measured by 
the observer Oq decrease, as it can be verified on Fig. EC! 



54. Effect on M, 



II) i.X OjTid - - 1 1 1 , ; v 



Table |2| gives (i) maximum mass rotating configurations 
without any magnetic field (first line), (ii) along a se- 
quence of constant magnetic dipole moment M (second 
line) and (iii) rotating configurations with magnetic fields 
close to the maximum value (third line), for the EOS Pol2, 
PandN and BJI. 

As in the static case, the maximum gravitational mass 
along sequences at fixed Ai, M max , increases with M. (cf. 
Figs. and E3|). A difference with the static case is that 




Fig. 22. Maximum Keplerian frequency Qk/2tt along rotating 
sequences at fixed magnetic dipole moment as a function of the 
magnetic dipole moment for the PandN EOS. 



20 30 
Mognetic moment [10 30 A.m 2 ] 



Fig. 24. Same as Fig. B2 but for the BJI EOS. 



the maximum baryon number is an increasing function of 
A4 for all EOS, including BJI (compare Figs. |^ and |23| ). 

For a given central-enthalpy and magnetic dipole mo- 
ment, the Keplerian angular velocity Qk is the value of 
f2 above which no stationary solution exists, the star be- 
ing at the break up limit under the effect of centrifugal 
inertial forces. The maximum Keplerian velocity along a 
sequence at fixed M. is reached at, or very close to, the 
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Table 2. Neutron star configurations in the rotating case. For each EOS, the first two lines are maximum mass configurations 
at fixed M. Q is the angular velocity; for the first two lines of each EOS £7 is the Keplerian velocity Ok- P = 2n/Sl is the 
rotation period, _R c irc the equatorial circumferential radius, J the angular momentum, Q the quadrupole moment as defined 
in Salgado et al. 1994a, U cq the fluid velocity at the equator as measured by a locally non-rotating observer and the 
"dragging of the inertial frames" at the star's centre. The remaining symbols are defined in the caption of Table We use 
G = 6.6726 1CT 11 m 3 kg -1 s -2 , c = 2.997925 10 s ms" 1 and 1 M© = 1.989 10 30 kg. 
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Fig. 21. Maximum gravitational mass (left) and maximum baryon number (right) along rotating sequences at fixed magnetic 
dipole moment as a function of the magnetic dipole moment for the PandN EOS. 
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Fig. 23. Same as Fig. |l| but for the BJI EOS. 

maximum mass configuration. As seen from Figs. ^2] and 
, the maximum SIk increases with M. 

6. Magnetized sequences at constant baryon num- 
ber 

Constant baryon number sequences may represent the 
time evolution of a neutron star and have been investi- 
gated in recent studies of non-magnetized rotating neu- 
tron stars (Cook et al. 1994b, Salgado et al. 1994a, b). 

We analyze in this section the behaviour of constant 
baryon number sequences at a fixed magnetic dipole mo- 
ment M. — 1.5 10 32 Am 2 , by stressing the differences with 
the non-magnetized case. 

6.1. Normal sequence 

A normal sequence is a sequence with a baryon num- 
ber lower than the maximum one supported by non- 
rotating configurations. We consider a sequence at bary- 
onic mass B = 3.00 M , built with the Pol2 EOS and the 
choice ( [33] ) for the current function /. Let us first men- 
tion that the choice M = 1.5 10 32 Am 2 corresponds to 
highly magnetized objects, the polar value of the mag- 
netic field varying from B = 1.4 10 3 GT (SI = 0) to 
B = 5.0 10 2 GT (SI = Sl K ) along the sequence. The dif- 
ference between various characteristic quantities of mag- 
netized and non-magnetized stars of the sequence are rep- 
resented in Fig. |2^. According to Fig. |25| a the magnetic 
field increases the gravitational mass, but only by a small 
amount (for the chosen value of Al), around one thou- 
sandth. The angular momentum increase is more substan- 
tial, being of the order of 3 % (Fig. f2ojb). The star's radius 



increases under the centrifugal effect of the Lorentz forces 
by roughly 5 % (Fig. ^B|c) and consequently the central en- 
ergy density decreases, by roughly 4 % (Fig. p5|c). For this 
normal sequence, the response of the global quantities to 
the magnetic field is the same as in the static case (cf 
Sect. |). 

We have also constructed normal sequences with the 
PandN EOS. Their behaviour is very similar to that de- 
scribed above. 

6.2. Supramassive sequence 

A supramassive sequence is a sequence with a baryon num- 
ber in excess of the maximum value that can be supported 
in the absence of rotation. We consider a supramassive 
sequence which in addition presents the phenomenon of 
spin-up by angular momentum loss (Cook et al. 1994a): 
the B = 3.80 M Q sequence built on the Pol2 EOS. The fact 
that along this sequence a decrease of angular momentum 
leads to an increase of angular velocity is demonstrated in 
Fig. 7 of Salgado et al. (1994a). The sequence begins at 
the angular velocity SI = 3.20 10 3 rads^ 1 and ends at the 
stability limit SI = 3.52 10 3 rads -1 . The electromagnetic 
effects are reported in Fig. ^6] and appear to differ from 
those of the normal sequence. The gravitational mass de- 
creases for the lowest angular velocities and increases for 
the highest ones (Fig. |26|a), whereas the angular momen- 
tum decreases all the time (Fig. p^b). This decrease in 
angular momentum can be attributed mainly to the star's 
radius decrease (Fig. |26|c). This latter behaviour contrasts 
with the normal sequence case (cf. Fig. p5|c) and can be 
explained as follows. Near the minimum angular velocity 
(SI = 3.20 10 3 rads -1 ), the centrifugal forces take their 
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0.5 1 1.5 

a [ 10 3 rad s" 1 ] 





0.5 1 1.5 

D [ 10 3 rad s" 1 ] 



Fig. 25. Differences between configurations without magnetic field and with a magnetic dipole moment M — 1.5 10 32 A m 2 along 
a sequence at constant baryonic mass B = 3.00 Mq for the Pol2 EOS and the current function f(x) = const: (a) gravitational 
mass difference AM; (b) angular momentum difference A J; (c) circumferential radius difference AR c i rc ; (d) central energy 
density difference Ae c . The path along the sequence is parametrized by the angular velocity fl. 



origin, not in a high rotation rate, but in a large spa- 
tial extension. The Lorentz forces play then a greater role 
in the equilibrium of the star since they are mainly lo- 
calised in the star's centre and are more efficient in sup- 
porting the star than the centrifugal forces, which are lo- 
calised at the star's periphery. The replacement of the 
centrifugal forces by the Lorentz ones in balancing grav- 
ity leads to a smaller stellar radius. This can be illus- 
trated by the following example: a star (built upon the 
Pol2 EOS) of baryon mass 3.77 M© and central density 
3.1 Pnuc has a circumferential radius of 26.5 km if it is 
supported by rotation (f2 = 3.3 x 10 3 rads _1 ) and only 
21.2 km if it is static and supported by the magnetic field 
(A4 = 7.02 x 10 32 Am 2 ), for the same number of baryons. 
This phenomenon explains the negative value of Ai? c ; rc 
in Fig. H= which reaches -9% at Q = 3.20 10 3 rads" 1 . 
When f2 increases the centrifugal forces become more effi- 
cient and the decrease of the star's radius with respect to 



the non-magnetized case becomes less pronounced. Due to 
the contraction of the star, the central density increases 
(Fig. p6|d). According to this explanation, the same phe- 
nomenon should exist in the Newtonian case, but we have 
not conducted a systematic study of Newtonian configu- 
rations. 



7. Conclusion 

We have extended an existing numerical code for comput- 
ing perfect fluid rotating neutron stars in general relativity 
(BGSM, Salgado et al. 1994a, b) to include the electromag- 
netic field. This latter is calculated by solving the relativis- 
ts Maxwell equations with an electric current distribution 
which is compatible with the star's equilibrium (i.e. the 
Lorentz force acting on the conducting fluid shall be the 
gradient of some scalar in order to balance gravity and 
the inertial centrifugal force). In order to preserve the sta- 
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a 




Q [ 10 3 rad s 1 ] 



c 




[ 10 3 rad s 1 ] 



Fig. 26. Same as Fig. [25| but for a supramassive sequence at I 

tionarity, axisymmetry and circularity properties of space- 
time, we consider only axisymmetric poloidal magnetic 
fields. The equations are numerically solved by means of 
a pseudo-spectral technique which results in a high accu- 
racy as tests on simple electromagnetic configurations (for 
which an analytical solution is available) have shown: the 
relative error on the electromagnetic field is of the order 
of 10~ 5 inside the star and 10 -9 outside it. The part of 
the code relative to the deformation of the star by Lorentz 
forces has been tested by comparison with Ferraro's ana- 
lytical solution in the Newtonian case. The fact that the 
numerical output is a solution of Einstein equations has 
been tested by two virial identities: GRV2 and GRV3. 

We have then used the code to investigate the effect of 
the magnetic field on rotating neutron stars. For this pur- 
pose we considered magnetic field amplitude ranging from 
zero up to huge values, of the order of 10 5 GT, which is 
ten thousand times bigger than the highest measured val- 
ues at the pulsars surfaces and is the value for which the 
magnetic pressure equals (with an opposite sign along the 
symmetry axis) the fluid pressure near the centre of the 



b 




fi [ 10 3 rad s 1 ] 



d 




Q [ 10 3 rad s 1 ] 



= 3.80 M Q . 

star. Let us note that such enormous magnetic fields are 
expected to decay on very short time scales via the mecan- 
ism of ambipolar diffusion, as investigated by Haensel et 
al. (1991), Goldreich & Reisenegger (1992) and Urpin & 
Shalybkov (1995), which is very efficient when the electric 
current is perpendicular to B, as in the present case. The 
decay time scale for non-superconducting matter com- 
puted by the above authors is 10 2 yr for B ~ 10 4 GT 
and 10 6 yr for B ~ 1 GT. 

According to our study, the influence of the magnetic 
field on the star's structure is mostly due to Lorentz forces 
and not to the gravitational field generated by the elec- 
tromagnetic stress-energy. This may be understood once 
it has been realized that a magnetic field of 10 5 GT has an 
energy density of 0.25 /9 n ucC 2 , whereas the matter density 
at the centre of neutron stars is between 1 and 10 /9 n uc- 
Although the electromagnetic energy is much lower than 
the fluid mass-energy, the deformation of the star can be 
as dramatic as that of Fig. || because of the anisotropic 
character of the magnetic pressure, just as the anisotropic 
centrifugal forces can highly deform the star in the rotat- 
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ing case though the kinetic energy is much lower than the 
fluid rest mass. 

In static and slowly rotating cases, Lorentz forces 
stretch the star away from the symmetry axis. The de- 
formation is appreciable only for B > 10 2 GT. In highly 
relativistic situations (supramassive sequences), the effect 
of the magnetic field is instead to reduce the star's equa- 
torial radius (at fixed baryon number). 

The maximum poloidal magnetic field supported by 
neutron stars has a polar value between 4 x 10 4 and 1.5 x 
10 5 GT depending upon the EOS and the rotation state 
of the star. Let us recall that the magnetic field at the 
star's centre is two to four times higher than at the poles. 

The impact of the magnetic field on the maximum 
mass of neutron stars is very limited for magnetic fields of 
the order of 1 GT, whereas it is important for the mag- 
netic fields near the maximum value (~ 10 5 GT): in the 
static case, M max is increased by 13% to 29% — depending 
upon the EOS — with respect to non-magnetized configu- 
rations. In fact, the magnetic field reveals to be more effi- 
cient in increasing M max than the rotation, except for the 
EOS HKP, where the maximum mass in rotation (with- 
out any magnetic field) is 21.2% higher than M max for 
stationary configurations (Salgado et al. 1994a), whereas 
the static magnetized M max is only 13.3% higher. For the 
PandN EOS, the M max increase by both mecanisms are 
similar (<~ 16%), whereas for the Pol2, BJI and Diaz II 
EOS, the static magnetized M max lies above the rotating 
non-magnetized M max . 

In the future, we plan to study the stability of the 
magnetized configurations presented in this article. Two 
types of instabilities may be expected to occur for high 
values of B: (i) a pure electromagnetic instability towards 
another electric current - magnetic field distribution (of 
lower energy) and (ii) a non-axisymmctric instability re- 
sulting in a triaxial stellar cquilibirum shape, which would 
be the magnetic analog of the transition from the MacLau- 
rin spheroids to the tri-axial Jacobi ellipsoids for high ro- 
tational velocities of Newtonian incompressible bodies. 
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